rm(list=ls())

library(cjoint)
library(tidyr)
library(tibble)
library(stargazer)
library(jtools)
library(dplyr)
library(haven)
library(readstata13)
library(ggplot2)
library(interactions)
library(ggpubr)
library(ggplotify)

#########################################################################
# REPLICATION FILES: MILITARIZATION AND PERCEPTIONS OF LAW ENFORCEMENT  #
# JOURNAL: BJPS                                                         #
# FIGURE 6: CONTINUOUS INTERACTIONS                                     #
#########################################################################

#### STEP 1: READ DATA                              ####

militconj <- read.dta13("MilitConjoint.dta")
baselines <- list()

#### STEP 2: DECLARE THEME                          ####

theme_bw1 <- function(base_size = 13, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = 11, colour = "black",  hjust = .5 , vjust=1),
      axis.text.y =       element_text(size = 11 , colour = "black", hjust = 0 , vjust=0.5), # changes position of X axis text
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=1.5),
      axis.title.x =      element_text(size =9),
      panel.background = element_rect(fill = "#f7f7f7"),
      strip.background = element_rect(colour="black", fill="white"),
      legend.position="bottom"
    )
}

#### STEP 3: INTERACTION PLOT, MILITARY PRESENCE    ####

# 1. EFFECTIVENESS
lm_effev <- lm(zeffective ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)
jev<- johnson_neyman(lm_effev, 
                     pred = uniform, 
                     modx = countevents, 
                     alpha = .05,
                     insig.color="grey",
                     sig.color="black",
                     title="Effectiveness",
                     mod.range=c(0,101))

jev <- jev$plot + xlab("Confrontations") + 
  ylab("Military uniform slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jev <- ggplot_build(jev)
jev$data[[8]]$alpha <- 0
jev <- ggplot_gtable(jev)
jev  <- as.ggplot(jev)

plot(jev)

# 2. CIVIL LIBERTIES
lm_libev <- lm(zlib ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)
jlib<- johnson_neyman(lm_libev, 
                      pred = uniform, 
                      modx = countevents, 
                      alpha = .05,
                      insig.color="grey",
                      sig.color="black",
                      title="Civil Liberties",
                      mod.range=c(0,101))

jlib <- jlib$plot + 
  xlab("Confrontations") + 
  ylab("Military uniform slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jlib <- ggplot_build(jlib)
jlib$data[[8]]$alpha <- 0
jlib <- ggplot_gtable(jlib)
jlib  <- as.ggplot(jlib)

plot(jlib)

# 3. CORRUPTION
lm_correv <- lm(zcorrupt ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)

jcorr<- johnson_neyman(lm_correv, 
                       pred = uniform, 
                       modx = countevents, 
                       alpha = .05,
                       insig.color="grey",
                       sig.color="black",
                       title="Corruption",
                       mod.range=c(0,101))

jcorr <- jcorr$plot + 
  xlab("Confrontations") + 
  ylab("Military uniform slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jcorr <- ggplot_build(jcorr)
jcorr$data[[8]]$alpha <- 0
jcorr <- ggplot_gtable(jcorr)
jcorr  <- as.ggplot(jcorr)

plot(jcorr)

# 4. NEIGHBORHOOD
lm_neighev <- lm(zneighb ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)
jneigh<- johnson_neyman(lm_neighev, 
                        pred = uniform, 
                        modx = countevents, 
                        alpha = .05,
                        insig.color="grey",
                        sig.color="black",
                        title="Neighborhood",
                        mod.range=c(0,101))

jneigh <- jneigh$plot + 
  xlab("Confrontations") + 
  ylab("Military uniform slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jneigh <- ggplot_build(jneigh)
jneigh$data[[8]]$alpha <- 0
jneigh <- ggplot_gtable(jneigh)
jneigh  <- as.ggplot(jneigh)

plot(jneigh)

#### STEP 4: INTERACTION PLOT, HOMICIDES            ####

#1. EFFECTIVENESS
lm_effhom <- lm(zeffective ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)
jeffhom2<- johnson_neyman(lm_effhom, 
                          pred = uniform, 
                          modx = thom, 
                          alpha = .05,
                          insig.color="grey",
                          sig.color="black",
                          title="",
                          mod.range=c(0,133))

jeffhom2 <- jeffhom2$plot + xlab("Average Homicide Rate") + 
  ylab("Military uniform slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jeffhom2 <- ggplot_build(jeffhom2)
jeffhom2$data[[8]]$alpha <- 0
jeffhom2$data[[9]]$alpha <- 0
jeffhom2 <- ggplot_gtable(jeffhom2)
jeffhom2  <- as.ggplot(jeffhom2)

plot(jeffhom2)

#2. CIVIL LIBERTIES
lm_libhom <- lm(zlib ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)
jlibhom2<- johnson_neyman(lm_libhom, 
                          pred = uniform, 
                          modx = thom, 
                          alpha = .05,
                          insig.color="grey",
                          sig.color="black",
                          title="",
                          mod.range=c(0,133))

jlibhom2 <- jlibhom2$plot + 
  xlab("Average Homicide Rate") + 
  ylab("Military uniform slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)

plot(jlibhom2)

#3. CORRUPTION
lm_corrhom <- lm(zcorrupt ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)

jcorrhom2<- johnson_neyman(lm_corrhom, 
                           pred = uniform, 
                           modx = thom, 
                           alpha = .05,
                           insig.color="grey",
                           sig.color="black",
                           title="",
                           mod.range=c(0,133))

jcorrhom2 <- jcorrhom2$plot + 
  xlab("Average Homicide Rate") + 
  ylab("Military uniform slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)

plot(jcorrhom2)

#4. NEIGHBORHOOD
lm_neighhom <- lm(zneighb ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)

jneighhom2<- johnson_neyman(lm_neighhom, 
                            pred = uniform, 
                            modx = thom, 
                            alpha = .05,
                            insig.color="grey",
                            sig.color="black",
                            title="",
                            mod.range=c(0,133))

jneighhom2 <- jneighhom2$plot + 
  xlab("Average Homicide Rate") + 
  ylab("Military uniform slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)

plot(jneighhom2)


#### STEP 5: GROUP INTERACTION PLOTS               ####

preshom <- ggarrange(jev, jlib, jcorr, jneigh, jeffhom2, jlibhom2, jcorrhom2, jneighhom2, ncol=4, nrow=2,
                       common.legend = TRUE, legend = "bottom")

plot(preshom)

#ggexport(preshom, 
         #filename="Figure6.png",
         #width = 900,
         #height = 600)

